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Abstract 


Currently,  the  Space  Defense  Center  (SDC)  uses  the  position  and 
velocity  data  of  a  satellite  track  to  update  the  orbital  elements  of 
the  satellite.  An  alternate  approach  would  be  to  process  the  radar 
data  at  the  tracking  site  and  transmit  two-body  orbital  elements  to  SDC 
as  measurement  updates.  This  would  significantly  reduce  the  data  load 
at  SDC.  The  truth  model  used  in  this  study  to  evaluate  filter  per¬ 
formance  and  provide  measurement  updates  to  the  filters  includes  the 
first  four  zonal  harmonics  and  the  first  sectoral  harmonic  of  the  geo¬ 
potential  and  lhe^air  drag.  The  estimator  models  include  two-body 

*  !  j  4-  •*  i 

dynamics  and  ^  perturbations.  In  addition,  the  derivative  of  the  semi¬ 
major  axis  is  used  to  estimate  the  effects  dueto  air  drag.  Perfect 
dynamics  with  measurement  noise  are  assumed  for  the  filter  models. 

Three  filters  are  evaluated;  Least  Squares,  Bayes,  and  Kalman.  The  per- 

i 

formance  of  the  three  filters  is  similar.  The  Bayes  and  Kalman  Filters 
are  nearly  the  same  and  indistinguishable  in  this  study.  All  three 
filters  performed  adequately  with  the  following  exceptions.  The  filters 
diverge  when  singularities  are  present  in  the  classical  orbital  elements. 
All  three  filters  underestimate  the  covariance  of  the  estimates. 

Finally,  the  filters  track  noise  in  the  estimate  of  the  derivative  of 
the  semi -major  axis,  degrading  the  prediction  capabilities  of  the 
filters. 
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ORBITAL  ESTIMATION  USING  TWO-BODY  CLASSICAL 


ORBITAL  ELEMENTS  AS  MEASUREMENT  UPDATES 

I .  Introduction 

The  Space  Defense  Center  (SDC)  located  in  the  NORAD  Cheyenne 
Mountain  Complex  (NCMC)  is  responsible  for  maintaining  current  orbital 
element  sets  for  all  earth  orbiting  satellites.  Currently,  there  are 
more  than  4500  satellites  orbiting  the  earth.  SDC  accomplishes  its 
mission  by  receiving  and  processing  data  from  various  spacetrack 
sensors  located  throughout  the  world.  The  data  is  processed  at  SDC  and 
the  satellite's  orbital  elements  are  estimated  at  SDC.  The  orbital 
elements  are  then  used  to  predict  the  satellite's  position  for  future 
acquisition  and  tracking. 

Currently,  SDC  estimates  the  orbital  elements  of  a  satellite  using 
the  satellite's  position  and  velocity  data  received  from  the  spacetrack 
sensors.  Either  a  least  squares,  batch  filter  or  a  Bayesian,  recursive 
filter.  Is  used  in  the  estimation  program.  An  average  satellite  orbits 
the  earth  between  twelve  and  sixteen  times  per  day.  Depending  upon  the 
satellite's  ground  track  and  SDC  tracking  requirements,  a  satellite  may 
be  tracked  10-15  times  dally.  Each  track  may  contain  fifteen  separate 
measurements  of  position  and  velocity,  or  ninety  individual  pieces  of 
data  per  track.  This  could  result  in  as  many  as  1000  measurements  per 
satellite  per  day  being  transmitted  from  the  sensors  to  SDC.  Also, 
this  data  must  be  stored  at  SDC  for  a  set  period  of  time  to  be 
available  for  estimating  the  satellite's  orbit  using  a  batch  filter. 


As  the  number  of  earth  satellites  increases  and  new,  improved  space- 
track  sensors  provide  more  and  better  data,  the  communications  lines 
to  SDC  and  the  data  storage  facilities  at  SDC  are  becoming  overloaded. 

The  purpose  of  this  thesis  Is  to  Investigate  the  possibility  of 
reducing  the  amount  of  data  transmitted  to  SDC  per  satellite  track  and 
the  amount  of  data  stored  at  SDC.  The  proposal  is  to  have  the  In¬ 
dividual  tracking  sites  process  the  position  and  velocity  on  each 
track  and  generate  two-body  orbital  elements  using  on-site  computers. 

The  site  will  then  transmit  the  two-body  orbital  elements,  six  pieces 
of  data  per  track,  to  SDC.  SDC  will  then  process  the  orbital  elements 
as  measurement  updates  to  an  estimator  which  would  include  pertubations 
to  the  two-body  orbit  in  the  estimate. 

The  approach  taken  in  this  thesis  is  divided  into  three  parts. 

First,  a  truth  model  will  be  developed  to  simulate  an  actual  satellite 
orbit  and  to  provide  data  for  the  estimators.  Next,  three  separate 
estimators  will  be  designed  to  estimate  the  satellite's  orbital  elements. 
Finally,  a  performance  analysis  of  the  three  estimators  will  be  con¬ 
ducted.  All  estimators  will  be  designed  to  accept  classical  orbital 
elements  as  measurements.  Various  satellite  orbits  will  be  investi¬ 
gated  to  evaluate  effects  due  to  atmospheric  drag,  high  and  low 
eccentricity,  low  inclination,  and  a  lack  of  tracking  data  due  to 
unusual  orbital  geometry. 

During  this  study  the  following  assumptions  are  made: 

1.  The  satellite  is  non-thrusting  and  unable  to  maneuver. 

2.  The  satellite  has  a  constant  drag  coefficient,  constant 
mass,  and  constant  surface  area. 

3.  The  geocentric  equatorial  coordinate  frame  is  considered  to 


'  HP 


be  an  Inertial  reference  frame. 

4.  The  earth's  gravitational  field  is  adequately  modeled  using 
zonal  harmonics  through  J4  and  the  first  sectoral  terms,  C22  and  S22* 

5.  The  satellite  dynamics  are  known  perfectly. 

6.  The  noise  in  the  sensor  tracking  data  Is  zero-mean,  white, 
Gaussian  noise. 

7.  The  atmospheric  density  above  at  an  altitude  greater  than 
800km  is  assumed  to  be  zero. 

This  thesis  is  presented  in  the  following  chapters.  Chapter  II 
describes  the  truth  model,  including  satellite  dynamics,  coordinate 
frame  used,  and  the  method  of  obtaining  measurements  for  use  in  the 
estimators.  The  third  chapter  presents  the  basic  estimator  model  and 
describes  the  least  squares,  Bayes,  and  Kalman  filters  used.  Chapter 
IV  describes  the  various  satellite  orbital  cases  used  for  filter  eval¬ 
uation.  The  fifth  chapter  presents  the  estimation  results  and  filter 
performance  evaluation.  The  sixth  chapter  contains  conclusions  of  this 
study  and  recommendations  for  future  work. 
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II .  Truth  Model 


The  truth  model  fulfills  two  Important  functions  In  this  study. 
First,  It  simulates  an  actual  satellite  orbit.  The  orbit  Is  simulated 
by  numerically  Integrating  the  equations  of  motion  for  the  satellite 
using  Cowell's  method.  In  order  to  do  this  the  coordinate  frame  and 
system  differential  equations  must  be  specified.  Before  integration, 
the  initial  conditions  for  the  differential  equations  must  also  be 
specified.  The  equations  of  motion  contain  the  two-body  accelerations 
plus  perturbative  accelerations  due  to  the  earth's  geopotential  and 
atmospheric  drag.  After  integrating  the  equations  of  motion,  the 
satellite's  position  and  velocity  are  output  from  the  truth  model. 
Zero-mean,  white,  gaussian  noise  is  then  added  to  the  position  and 
velocity  data,  simulating  actual  radar  tracks.  Finally,  the  position 
and  velocity  data  are  converted  to  the  classical  orbital  elements  for 
use  as  measurement  updates  to  the  estimators.  A  description  of  the 
orbital  simulation  and  measurement  updates  will  be  completed  in  the 
following  paragraphs. 

Coordinate  Frame 

Before  integrating  the  equations  of  motion  for  a  satellite,  a 
suitable  inertial  reference  frame  must  be  defined.  Additional  reference 
frames,  some  non-inertial ,  may  be  defined  to  aid  in  computation. 

Three  coordinate  frames  are  used  in  this  study;  geocentric-equatorial, 
perifocal,  and  local  horizon  coordinate  frames. 

For  earth  orbiting  satellites,  the  geocentric  equatorial  frame  can 
be  assumed  to  be  an  inertial  reference  frame.  As  shown  in  Figure  1, 
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the  center  of  this  frame  is  the  center  of  the  earth  and  the  principle 
plane  is  the  earth's  equatorial  plane.  The  positive  X-axis  points  in 
the  direction  of  the  vernal  equinox  and  the  positive  Z-axis  points  in 
the  direction  of  the  north  pole.  Finally,  the  Y-axis  completes  the 
right-handed  XYZ  coordinate  frame.  It  can  be  seen  in  Figure  1  that 
the  XYZ  system  is  not  rotating  with  the  earth,  but  is  fixed  with 
respect  to  the  stars.  This  system  is  used  primarily  for  expressing 
and  integrating  the  equations  of  motion. 

The  perifocal  coordinate  frame  is  shown  In  Figure  2.  The  funda¬ 
mental  plane  for  this  system  is  defined  by  the  plane  of  the  satellite's 
orbit.  The  principle  axis,  Xp,  points  toward  the  perigee  of  the  orbit. 
The  Yp  axis  is  rotated  90°  in  the  direction  of  orbital  motion  in  the 
orbital  plane.  The  Zp  axis  completes  the  right-handed  XpYpZp  system. 
The  perifocal  system  is  used  for  transforming  from  classical  orbital 
elements  to  satellite  position  and  velocity. 

The  final  coordinate  frame  is  the  local  horizon  coordinate  frame. 
The  center  of  this  frame  is  at  the  intersection  of  the  satellite's 
position  vector  from  the  center  of  the  earth  and  the  surface  of  the 
earth.  The  directions  of  this  system  are  up  along  the  radius  vector. 
East,  and  North.  This  frame  is  convenient  for  expressing  the  earth's 
geopotential . 

In  the  course  of  this  study.  It  will  be  necessary  to  transform 
from  one  coordinate  frame  to  another.  The  direction  cosine  matrices 
necessary  for  the  various  transformations  are  shown  in  Appendix  A. 
States  -  State  Equations 

The  system  states  used  in  the  truth  model  are  the  satellite 
position  and  velocity  in  the  geocentric-equatorial  reference  frame. 


vernal  equinox  direction 
Figure  1.  Geocentric  Equitorial  Frame 


Figure  2.  Perifocal  Coordinate  Frame 


The  state  equations  are  shown  In  equation  (1), 


(1) 

/ 


where  ax.  ay,  and  a2  are  accelerations  due  to  perturbative  forces. 

These  forces  will  be  discussed  in  later  sections. 

The  units  used  In  these  equations  are  geocentric,  canonical  units. 
In  this  system  the  distance  unit,  DU,  Is  the  earth's  mean  equatorial 
radius,  equation  (2). 

1  DU  =  6378.145  Km  (2) 

The  time  unit,  TU,  is  chosen  so  that  y  is  equal  to  one. 

1  TU  =  13.44686457  min  (3) 

Initial  Conditions 

Initial  conditions  are  required  for  the  satellite's  position  and 
velocity  in  order  to  integrate  the  equations  of  motion.  These  condi¬ 
tions  are  obtained  by  specifying  a  particular  orbit  at  a  particular 
time.  The  orbital  elements  are  transformed  to  a  position  and  velocity 
vector  In  the  perifocal  coordinate  frame.  The  position  and  velocity 
vector  are  then  transformed  to  the  geocentric-equatorial  coordinate 
frame,  providing  initial  conditions  for  the  state  equations  given  In 
equation  (1).  The  technique  for  obtaining  the  position  and  velocity 
from  classical  orbital  elements  is  shown  in  Appendix  B. 
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Earth's  Geopotential 

The  gravitational  potential  of  the  earth  at  a  point  at  a  distance 
r  from  the  center  of  the  earth  is, 

U  =  J  (1  +  1  (^)  ?  p£  (sin  6)  [Cnm  cos  mx  +  Snm  sin  mx]  (4) 
n=2  m=2 

where, 

p  is  the  earth's  gravitational  constant 

re  is  the  earth's  mean  equatorial  radius 

p”  is  the  associated  legndre  function  of  degree  n  and  order  m 
Cnm  &  Snm  are  numerical  coefficients 

6  is  the  geocentric  latitude  of  the  point 

X  Is  the  geocentric  longitude  of  the  point 

For  m  =  0,  the  Snm  coefficients  are  not  needed  and  the  Cnm  coef¬ 

ficients  are  replaced  by 


Jn  =  -Cno  (5) 

These  harmonics  are  independent  of  longitude  and  are  called  zonal 
harmonics.  The  zonal  harmonic,  J] ,  is  set  equal  to  zero  by  choosing  a 
coordinate  frame  with  the  origin  at  the  earth's  center  of  mass.  For 
this  study,  only  the  terms  through  J4  were  retained.  The  terms  for 
m  =  n  are  called  sectoral  harmonics.  Only  the  first  sectoral  harmonic, 
m  =  n  =  2,  Is  used  In  this  study.  The  remaining  terms,  called  tesseral 
harmonics,  are  not  Included  In  this  study. 

The  first  term  In  equation  (4)  Is  the  major  term  In  the  equations 
of  motion  for  a  satellite  and  represents  the  two-body  motion  of  the 
satellite.  This  term  Is  represented  In  equation  (1)  by  the  first 
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acceleration  terms.  The  remaining  terms  In  the  geopotential  are 
Included  in  the  perturbation  term. 

The  terms  of  the  acceleration  due  to  the  earth's  geopotential  are 
most  easily  expressed  In  the  local  horizon  coordinate  frame  [Ref  9, 
88-92].  In  this  frame 


M  A  A  A 

'^3  +  a  =  9u  u  +  9e  e  +  9n  n 


(6) 


where 


„  _  3U 
gu  ■  af 


-  +  z  (n+1 )  £  "  z  P*  (sin5)[C 

r2C  n=2  r  m=0 


.nm  cos  m  +  Snm  sin  mx]j" 
(7) 


9e 


1  3U 


r  cos6  9X 


=  -u  |  £e  n  "  m  P™  (s1n6)[Cnm  sin  mx  -  Snm  cos  mx] 

r2  cos6  n=2  r  m=0  ^ 

a  -  1  9U 

9n  ■  F  36 

=  _jj  ?  !ja  "  COS6  3  (P„  (sine))  [Cnm  cos  mx+Snm  sin  mx] 

~  r2-n-o  r  *  36 

n-2  m-0 


The  acceleration  terms  are  then  transformed  to  the  geocentric- 
equatorial  frame  for  integration.  The  zonal  and  sectoral  harmonics 
and  the  associated  legendre  polynomials  required  for  this  study  are 
presented  In  Appendix  C. 

Atmospheric  Drag 

The  acceleration  on  a  satellite  caused  by  atmospheric  drag  is 


a 


Bpva 


(10) 
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where. 


B  =  Cn  —  Is  the  satellite's  ballistic  coefficient 
u  m 

Co  ■  Is  the  drag  coefficient  of  the  satellite 
A  Is  the  satellite's  cross-sectional  area  perpendicular  to  the 
velocity 

m  is  the  satellite's  mass 

p  is  the  atmospheric  density  at  the  altitude  of  the  satellite 

va  Is  the  velocity  of  the  satellite  relative  to  the  atmosphere. 

The  satellite's  ballistic  coefficient  is  assumed  to  be  constant 
for  this  study  and  is  input  into  the  system  with  the  initial  orbital 
elements . 

The  velocity  of  the  satellite  relative  to  the  atmosphere  is 


x  +  o>y 


z 


(ID 


where 

x,  y,  and  z  are  the  satellite's  geocentric-equatorial  coordinates 
ui  =  .0588336001  -=^j-  is  the  earth's  angular  rotation. 

The  atmospheric  density,  p  ,  must  be  calculated  from  a  model 
atmosphere.  The  model  atmosphere  in  this  study  was  developed  by 
Jacchia  in  1960  [Ref  7].  The  model  accounts  for  altitude  above  a 
non-spherical  earth  and  diurnal  variations  In  the  density.  The  density 
is  given  by 

hO  Cl  +  0.19  [exp(0.0055H)  -  1.0]  [L+  cos.lj3(  (12) 

p  "  p0  100  C  £  s  m3 
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where 


Log  10  P0  =  -16.021  -  0.001985H  +  6.363  exp(-0.0026H) 


(13) 


H  =  6378.145  [r-1  +  f  u\  +  f  f2  u*  (1-uf)]  (14) 

is  the  height  of  the  satellite  above  the  earth, 
r  is  the  satellite's  radius 

f  =  1/298.3  is  the  earth's  flattening 
u2  =  sin  6 

6  is  the  satellite's  geocentric  declination 

Fio  is  the  solar  flux  at  10.7  cm 

cosiji  =  ux  cos  (0-50°)  +  Uy  sin  (0-50°) 

D  =  (Day  of  year)  X  360°/365.25  days 

A  A  A  A  A 

ux=  u  •  1  /  u  Is  a  perifocal  frame  unit  vector,  i  and  j  are 

A  A  7 

Uy  -  u  •  j  \  geocentric  equatorial  unit  vectors. 

For  this  study  it  is  assumed  that  p  =  0  for  altitudes  greater  than  800km. 
Measurements 

After  the  equations  of  motion  are  integrated,  orbital  elements 
must  be  calculated  for  updates  to  the  estimators.  The  satellite's 
geocentric-equatorial  position  and  velocity  are  output  from  the 
numerical  integration.  Zero-mean,  white,  gaussian  noise  is  then  added 
to  the  data  to  simulate  errors  in  the  radar  data.  The  standard 
deviation  is  500m  in  the  position  and  lOm/sec  in  the  velocity.  The 
position  and  velocity  data,  including  noise  is  then  converted  to 
classical  orbital  elements.  These  elements  are: 

a  -  semi -major  axis 
e  -  eccentricity 
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.1  |l  I  ||  1  .  JHppjS 


1  -  Inclination 

a)  -  argument  of  perigee 

Q  -  right  ascension  of  ascending  node 

M  -  mean  anomaly 

t  -  time  of  track 

The  method  of  transforming  from  a  position  and  velocity  vector  to 
classical  orbital  elements  is  shown  In  Appendix  D. 

i 
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III.  Estimation 


In  orbit  determination,  the  system  model  Is  a  non-linear  set  of 
equations.  The  parameters  in  the  model  are  well  known  and  the  assump¬ 
tion  of  perfect  dynamics  is  adequate  for  the  estimation  problem.  Also 
the  observation  relations  are  in  general  non-linear.  The  observations 
are  corrupted  by  zero-mean,  white  noise. 

Three  separate  estimators  are  used  to  estimate  the  orbital  param¬ 
eters.  The  filters  are  a  least  squares  filter,  a  Kalman  filter,  and  an 
inverse  covariance  Bayes  filter. 

The  model  dynamics,  observation  relationship,  and  the  three  filters 
are  discussed  in  the  following  sections. 

Dynamics 

The  dynamics  In  this  section  are  the  same  as  the  dynamics  discussed 
in  Chapter  II;  however,  a  simpler  model  is  proposed  to  reduce  computa¬ 
tional  time  and  expense.  The  model  proposed  in  this  section  is  essen¬ 
tially  the  same  as  the  model  currently  used  at  SDC  in  their  standard 
model  [Ref  14]. 

The  system  models  of  the  estimators  are  designed  to  model  the  two- 
body  dynamics  plus  effects  and  nlr  drag.  Only  the  secular  terms  in 
O2  are  considered.  The  secular  effects  due  to  J2  appear  in  the  deriv¬ 
atives  of  argument  of  perigee,  right  ascension,  and  mean  anomaly.  The 
presence  of  atmospheric  drag  reduces  the  semi -major  .axis  a  and  eccen¬ 
tricity  e.  Due  to  the  difficulty  in  modeling  the  atmosphere,  and  there¬ 
fore  the  air  drag,  the  derivatives  of  a  and  e  are  Included  In  the  system 
states  and  the  atmospheric  effect  is  estimated.  The  system  states  used 
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In  the  estimators  are  (see  List  of  Symbols): 


xT  =  [a,  e,  1,  a,  n,  M,  a,  e]  , 


(15) 


An  approximate  solution  to  the  non-linear  system  Is  given  by 


or. 


where 


x(t)  =  f  (x  (to),  t) 

(16) 

a(t)  =  a0  +  a ( t-t0) 

(17) 

e(t)  =  e0  +  e(t-t0) 

(18) 

i(t)  =  i0 

(19) 

Ol(t)  =  (DO  +  dt(t-t0) 

(20) 

n(t)  =  n0  +  n(t-t0) 

(21) 

M(t)  =  M0  +  n|c(t-t0)  +  \  n(t-t0)2 

(22) 

•  • 
a  =  ao 

(23) 

e  »  e0 

(24) 

a  =  <2  si"z  1  '  2> 

2a20-e2)2  2 

(25) 

a  -  -  cos  i 

2a2(l-e2)2 

(26) 

#  =  f  -5/2. 

2a  a 

(27) 

n|<  =  n  +  e 

(28) 

is  the  Kozai  mean  motion 

;  =  d  sin*  1-1) 

2a2(1.e2)3/2  2 

(29) 
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A  linearization  of  equations  (17)  through  (24)  yields  the  state 


transition  matrix,  <j>,  an  8  X  8  matrix.  4>  Is  given  by 

4>  =  V'f 


(30) 


The  elements  of  $  are  given  In  Appendix  E. 

Observation  Relationships 

The  measurement  updates  to  a  filter  are  given  in  general  by 


z  =  h  (x(t0) »  t)  +  v 


(31) 


where  h^  must  be  determined  through  geometric  relationships  between  the 
measurements  and  the  states.  For  example,  if  the  satellite's  position 
and  velocity  vectors  are  used  as  measurement  updates,  then  a  set  of  v 
equations  similar  to  those  in  Appendix  B  would  be  required.  Also 
coordinate  transformations  from  the  station  centered  to  perifocal  ref¬ 
erence  frames  would  be  required.  This  is  currently  the  method  used  at 
SDC .  In  this  study  the  measurements  are  the  classical  orbital  elements. 
These  measurements  are  a  linear  combination  of  the  states. 


z  =  H  x  +  v 


where 


v  is  a  zero-mean,  white  noise. 


(32) 


(33) 
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V 


The  associated  covariance  matrix  for  the  noise  v.  Is  Q.  It  Is 
assumed  for  this  study  that  the  noise  components  are  uncorrelated,  so 
that  the  off-diagonal  components  of  Q  are  zero. 

In  order  to  determine  the  diagonal  elements  of  Q,  it  is  assumed 
that  the  initial  tracking  data  had  standard  deviations  of  500m  in  the 
position  vector  and  lOm/sec  in  the  velocity  vector.  This  noise  was 
added  to  the  position  and  velocity  vector  and  the  orbital  elements  were 
calculated.  This  process  was  repeated  ten  times  and  then  an  average 
error  for  each  orbital  element  was  calculated.  The  results  show  that 
a  single  value  for  the  diagonal  elements  of  Q  is  adequate.  The  value 
used  in  this  study  is  10®. 

Least  Squares  Fil ter 

The  least  squares  filter  used  in  this  study  is  a  batch  filter. 

Thus,  every  time  an  update  is  required  the  data  must  be  batch  processed 
to  update  the  state  estimate.  In  this  study,  the  state  vector  is  up¬ 
dated  after  every  ten  measurements,  using  the  previous  estimate  as  the 
reference  trajectory.  The  basic  algorithm  used  in  this  study  is  given 
In  the  following  steps  [Ref  13,  62-63]. 

1.  Input  measurements,  zi  at  times  ti ,  measurement  noise  Qi ,  and 
reference  orbit  Xr(to).  zj  will  be  a  6  X  10  matrix  for  the  six  elements 
measured  at  ten  different  times.  As  stated,  Q  is  a  diagonal  matrix 
with  a  single,  constant  value  at  each  diagonal  element.  Also,  the 
reference  trajectory  Xr(to)  the  previous  state  estimate  evaluated 

at  time  to. 

2.  Propagate  the  state  vector  Xp  to  each  time  t]  and  calculate 
xjt-j)  and  the  state  transition  matrix  i>(ti,  t0). 
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3.  Assemble,  Q  and  H. 


Qi 


Q2 


Q10 


(34) 


(35) 


Hi  $(ti,  t0) 

H2  «(t2.  t0) 

H10  *(*10*  to^j 

4.  Calculate  the  covariance  matrix,  P. 

P  =  (H*t  Q-1  H*)-1  (36) 

5.  Calculate  the  residuals  r(ti)  and  assemble  into  a  vector. 

r(ti)  =  z(tj)  -  H(tj )  x(ti)  (37) 


r  = 


n(ti)~ 

l(t2) 

« 

« 

n(tio) 


(38) 


6.  Calculate  the  change  in  the  state  vector  <5x(t0). 

6x(t0)  =  P  H*t  Q'1  r 

7.  Calculate  the  estimated  state  vector  x_(to). 


(39) 


x(t0)  =  2ir(to)  +  ®2L(to) 


(40) 
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8.  Check  for  convergence.  For  this  study  6xi  is  compared  to 
*/FiT.  If  the  absolute  value  of  the  change  in  each  element  of  the  state 
vector  is  less  than  its  associated  standard  deviation,  the  solution  is 
assumed  to  have  converged  and  the  associated  estimate  is  assumed  to  be 
the  state  estimate  at  the  epoch  time,  to* 

9.  If  the  solution  has  not  converged,  the  new  solution  replaces 
the  reference  orbit.  This  new  reference  orbit  is  propagated  to  each 
time,  ti,  and  new  residuals  are  calculated  as  shown  in  Step  5.  This 
process  is  repeated  until  convergence  is  achieved  or  a  maximum  number 
of  iterations  is  exceeded.  The  maximum  number  of  iterations  used  in 
this  study  is  twenty. 

The  epoch  time,  t0,  in  this  filter  is  always  taken  as  the  time  of 
the  last  measurement. 

Since  this  filter  uses  only  the  last  ten  measurements,  it  is  a 
finite  memory  filter  and  should  not  encounter  problems  with  the  co- 
variance  matrix,  P,  becoming  singular.  A  larger  batch  size  or  more 
frequent  updates  to  the  estimates  may  be  implemented;  however,  a  standard 
batch  size  of  ten  measurements  per  update  is  maintained  in  this  study. 

Some  observations  concerning  the  filter  can  be  made.  First,  when 
using  only  the  last  ten  measurements,  the  dimension  of  some  of  the 
matrices  becomes  rather  large.  The  Q  matrix  is  a  60  X  60  matrix  in 
general;  however,  it  is  at  worst  block  diagonal  and  in  this  case 
strictly  a  diagonal  matrix.  The  P  matrix  is  an  8  X.8  matrix  which  must 
be  inverted.  Both  the  P  matrix  and  the  H*  matrix  are  calculated  once 
per  state  estimate  update  and  are  then  assumed  to  be  constant  for  that 
update  cycle.  A  more  exact  method  would  be  to  update  the  two  matrices 
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each  iteration.  For  this  problem  the  possible  decrease  in  the  number 
of  iterations  was  not  worth  the  corresponding  increase  in  computational 
time. 

Bayes  Fi 1  ter 

The  Bayes  Filter,  or  inverse  covariance  filter,  used  in  this  study 
is  a  recursive  filter,  providing  a  state  estimate  update  after  every 
measurement  input  into  the  filter.  This  filter  propagates  the  estimated 
state  vector  and  its  associated  covariance  forward  to  the  time  of  the 
current  measurement.  This  estimate  and  the  measurement  are  then  com¬ 
bined  to  form  a  new  estimate  of  the  state  vector  along  with  its  as¬ 
sociated  covariance  matrix.  The  algorithm  used  in  this  study  is  listed 
below  [Ref  13,  82-85]. 

1.  Input  the  previous  estimate,  x(ti  -  1),  new  data,  z_(t0)»  and 
its  noise  matrix,  Q. 

2.  Propagate  the  state  estimate  and  the  covariance  matrix  to  the 
time  of  the  new  measurement,  t^ . 

x(ti“)  =  f(x(ti-i),  tj  -  ti-i)  (41) 

P(ti')  =  $(ti , t-j _•] )  P(tj-i)  *T(ti ,  t^T )  (42) 

The  initial  state  vector  in  the  filter  is  the  first  measurement. 

The  initial  covariance  matrix  input  into  the  system  is 

Pij  =  1.  X  10"4  i  *  j  (43) 

Pii  =  1.  X  10'2  (44) 

3.  Calculate  the  new  covariance  matrix,  P(ti). 

P(ti )  =  {p-1 (ti“)  +  $T  HT  Q"1  H*}"1  (45) 
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4.  Calculate  the  gain  matrix,  G. 


G  =  P(ti)  $T  HT  Q"1  (46) 

5.  Calculate  the  residual  vector,  r. 

r  ■  z(ti)  -  Hx.( tn- ' )  (47) 

6.  Calculate  the  change  in  the  state  vector,  Sx. 

6x(t)  =  <5x(“)  +  G(r  -  H$Sx(-))  (48) 

starting  with 

Sx(-)  =  0  (49) 

7.  Check  for  convergence.  For  this  study  convergence  is  achieved 
if 


£  =  H$Sx(-)  <£  (50) 

where 

ei  =  1.  x  10'3  (51) 

If  the  solution  converges,  the  state  estimate  is 

x(t0)  =  x(ti~)  +  6x(t)  (52) 

8.  If  the  solution  has  not  converged,  new  residuals  are  calculated 
based  on  the  updated  state  estimate  and  steps  6  through  8  are  repeated. 
If  the  solution  has  not  converged  after  ten  iterations  the  algorithm 
is  exited  and  the  new  state  estimate  is  calculated  as  shown  in 
equation  (52). 
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The  data  storage  requirements  and  the  dimensions  of  the  matrices 
are  greatly  reduced  for  this  filter  when  compared  to  the  least  squares 
filter.  The  largest  dimension  of  a  matrix  in  this  filter  is  8X8. 

This  filter  provides  a  state  estimate  update  after  every  measurement. 
This  implies  that  the  state  estimate  is  more  current  than  the  least 
squares  filter;  however,  this  filter  must  generate  an  estimate  ten  times 
for  every  estimate  generated  by  the  least  squares  filter.  Like  the 
least  squares  filter,  the  gain  and  covariance  matrix  in  this  filter 
are  calculated  once  per  update  cycle.  These  matrixes  are  then  assumed 
to  be  constant  for  the  remaining  iterations.  Also,  like  the  least 
squares  filter,  the  filter  requires  one  8X8  matrix  inversion  per 
update  cycle. 

Kalman  Filter 

The  third  filter  investigated  on  this  study  is  a  Kalman  Filter. 

The  operation  of  the  Kalman  Filter  is  quite  similar  to  the  Bayes  Filter. 
It  is  a  recursive  filter  that  provides  a  new  state  estimate  following 
each  measurement  update.  The  algorithm  used  in  this  study  is  listed 
below  [Ref  13,  90-91]. 

1  and  2  -  same  as  Bayes  Filter. 

3.  Calculate  the  Kalman  Filter  gain,  K. 

K  =  P(tj-)  $T  Ht(Q  +  H$P(tj-}$  T  H1)-1  (53) 

4.  Calculate  the  new  covariance  matrix  P(t-j+) , 

P(ti+)  =  P(ti-)  -  KHsP(ti-)  (54) 
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5.  Calculate  the  residuals,  r. 

r  =  z  -  Hx(t-j-)  (55) 

6.  Calculate  the  change  in  the  state  vector,  6x(+). 

6x(+)  =  <5x(-)  +  K(r  -  H'jxSx(-))  (56) 

7.  Check  for  convergence.  The  convergence  check  is 

r  -  H$6x(")  <£  (57) 

where 

ei=1.0x  10“3  (58) 

If  convergence  is  achieved,  the  new  estimate  is 

x(ti+)  =  x(ti-)  +  6x(+)  (59) 


8.  If  convergence  is  not  attained,  calculate  new  residuals  using 
the  updated  state  vector  estimate  and  repeat  Steps  6  through  8.  Again, 
the  maximum  number  of  iterations  is  10. 

The  characteristics  of  the  Kalman  Filter  are  the  same  as  the  Bayes 
Filter  except  the  Kalman  Filter  requires  a  6  X  6  matrix  inversion 
instead  of  an  8  X  8  matrix  inversion. 

A  summary  of  the  three  filter  characteristics  is  listed  in  Table  I. 
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TABLE  I 

F i  1  tejr  Characteristics 

Characteristics 
Method  of  Update 
Frequency  of  Updates 

Storage  Requirements 

Dimension  of  Matrix 
Inversion 


Least  Squares 


Bayes 


Kalman 


Batch 

1  per  10  measure¬ 
ments 

Current  estimate 
and  covariance 
plus  last  10 
measurements 


Recursive 


Recursive 


every  measure-  every  measure¬ 
ment  ment 

Current  estimate  and  covar¬ 
iance  plus  last  measurements 


8X8 


8X8 


8X8 


IV.  Satellite  Cases 


Various  satellite  cases  were  investigated  in  this  study  to  eval<- 
uate  the  performance  of  the  three  filters.  The  orbital  elements  for 
each  case  are  listed  below  along  with  an  explanation  of  the  objectives 
of  each  case  (see  List  of  Symbols). 

Case  1 . 

a0  =  1.2  DU  =  7653.774  km 
e0  =  .1 

i0  =  .5  rad  =  28.648° 
w0  =  1.0  rad  =  57.296° 
n0  =  2.0  rad  =  114.592° 

M0  =  3.0  rad  =  171.887° 

At  =  10  TU  =  134.686  min 

where 

At  is  the  time  between  measurement  updates. 

The  satellite  orbit  is  used  to  demonstrate  the  operation  of  the 
filters  under  normal  operation  with  no  singularities  in  the  orbital 
elements.  Normal  operation  is  assumed  to  be  one  radar  track  every  two 

hours.  i  • 

j  / 

i 

Case  2. 

a0  =  2.0  DU  =  12,756.286  km 
e0  =  .1 

i0  =  1.0  rad  =  57.286° 
u)0  =  2.0  rad  =  114.592° 
fl0  =  3.0  rad  =  171.887° 
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M0  =  1.0  rad  =  57.296° 

At  =  50  TU  =  672.343  min 

The  objective  of  this  orbit  is  to  investigate  the  effect  on  the 
estimators  with  a  long  time  delay  between  updates.  There  are  approxi¬ 
mately  two  updates  per  day  with  this  orbit. 

Case  3. 


a0  =  1.15  DU  =  7335  km 

eo  =  .1 

i0  =  .5  rad  =  28.648° 
o)0  =  1.0  rad  =  57.296° 

=  2.0  rad  =  114.592° 

M0  =  3.0  rad  =  57.296° 

At  =  10  TU  =  134.686  min 

This  orbit  is  a  low  altitude  orbit.  This  orbit  should  identify 
problems  in  estimating  air  drag.  The  perigee  altitude  for  this  orbit  is 
200  km.  Also  the  drag  coefficient  of  the  satellite  was  increased  by 
a  factor  of  10,  from  0.01  to  0.1,  to  increase  the  effect  of  air  drag. 

Case  4. 

a0  =  1.2  DU  =  7653.774  km  | 

e0  =  0. 

i  o  =  0. 

U)q  =  0 . 

fi0  =  0. 

M0  =  0. 

At  =  10  TU  =  134.686  min 


25 


This  orbit  introduces  singularities  when  estimating  the  classical 
elements  from  position  and  velocity  data. 

Case  5. 

a0  =  4  DU  =  25500  km 
®o  •  ^ 

i0  =  1.1  rad  =  63° 
u)Q  =  4.7  rad  =  270° 

=  1-0  rad  =  57° 

M0  =  0  rad  =  0° 

At  =  10  TU  =  134.686  min 

This  orbit  simulates  a  highly  eccentric  orbit.  To  simulate 
tracking  problems  and  observation  problems,  all  data  generated  at 
altitudes  greater  than  6000  km  were  disregarded.  In  this  case  the  At 
given  will  vary,  depending  upon  the  altitude. 

The  performance  of  the  filters -using  each  of  the  orbits  is  pre¬ 
sented  in  Chapter  V.  The  filters  are  evaluated  for  twenty  measurement 
update  cycles,  which  is  two  days  for  the  high  frequency  tracks  and  ten 
days  for  the  low  frequency  tracks. 
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V.  Filter  Performance  Evaluation 


The  three  filters  used  in  this  study  are  evaluated  based  on  each 
filter's  estimate  of  the  covariance,  the  mean  error  of  the  state 
estimates,  and  the  root  mean  square  errors  of  the  state  estimates. 
Covariance  Analysis 

Since  the  dynamics  in  this  problem  are  non-linear,  a  Monte  Carlo 
analysis  is  required  to  evaluate  the  covariance  estimate  of  the  three 
filters.  The  first  step  is  to  determine  the  number  of  runs  required 
for  the  covariance  from  the  Monte  Carlo  analysis  to  reach  a  constant 
value.  The  variance  of  the  filter  is  given  by  the  following  equation: 


var(x) 


(60) 


/ 


where 

02  «  1  5  (XT  -  X/)2  (61) 

a  fTT  j=l  J 

N 

e  =  J  ^(XT-Xj)  (62) 

O 

ac  is  the  mean  squared  error 
e  is  the  mean  error 

N  is  the  number  of  runs  starting  from  2 
Xj  is  the  true  value  of  the  state 

A 

Xj  is  the  estimated  value  of  the  state 


This  analysis  provides  the  diagonal  elements  of  the  filter  covar¬ 
iance;  the  variance  of  the  state  estimates. 
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The  Bayes  and  Kalman  Filters  are  evaluated  using  twenty-five 
separate  runs  of  Case  I  of  the  satellite  orbits.  The  Least  Squares 
Filter  is  evaluated  using  twenty-five  runs  of  Case  II  of  the  satellite 
orbits.  The  mean  squared  error  is  calculated  at  the  time  of  the  last 
measurement  update. 

At  =  200  TU  *  45  hrs 

O 

The  mean  squared  error  a  for  each  state  is  plotted  as  a  function 
of  the  number  of  runs.  Figure  3  represents  a  plot  of  the  mean  squared 
error  of  the  inclination  versus  the  number  of  runs  for  the  Least 
Squares  Filter.  This  plot  shows  the  desired  characteristic  of  approach¬ 
ing  a  constant  value  after  about  12-14  runs.  Plots  for  the  semi -major 
axis  and  mean  anomaly  for  each  filter  are  shown  in  Figures  4-9.  The 
plots  for  the  remaining  states  are  included  in  Appendix  F.  As  can  be 
seen  in  the  plots  (i.e.,  Fig  6)  not  all  curves  look  like  Figure  3.  In 
some  cases  it  is  difficult  to  tell  when  and  if  the  mean  squared  error 
has  reached  a  nearly  constant  value.  However,  in  most  cases,  nearly 
constant  values  are  achieved  after  fifteen  runs.  Fifteen  runs  are 
used  for  all  subsequent  Monte  Carlo  analyses. 

Next,  the  variance  obtained  from  the  Monte  Carlo  analysis  is  com¬ 
pared  to  the  estimate  of  the  variance  from  the  filters.  The  variances 
of  the  semi -major  axis  and  mean  anomaly  are  used  for  comparison.  Table 
II  shows  a  comparison  of  the  variances  of  these  two  elements  for  each 
filter.  The  final  estimate  for  Case  I  is  used  for  comparison. 

As  shown  in  the  table,  all  three  filters  underestimate  the  variance 
of  both  elements.  Similar  results  are  obtained  for  the  other  four 
elements.  This  is  probably  due  to  the  assumption  of  perfect  dynamics 
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Fig  3.  o\  vs  N,  Least  Squares  Filter 


Fig  4.  oa6  vs  N,  Least  Squares  Filter 


Fig  5.  Oft  vs  N,  Least  Squares  Filter 


Fig  6.  oa  vs  N,  Bayes  Filter 


t.0I*  SIXB-A 


Fig  8.  oa2  vs  N,  Kalman  Filter 
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TABLE  II 


Covariance  Analysis 

Semi -major  Axis  (DU^)  Mean  Anomaly  (Rad2) 


Filter 

Estimate 

Monte  Carlo 

Estimate 

Monte  Carlo 

Least  Squares 

1.7  X  10"11 

6.3  X  lO"6 

4.7  X  lO’9 

1.4  X  lO"3 

Bayes 

2.2  X  10“12 

1.2  X  10'6 

2.6  X  lO’9 

5.3  X  10"5 

Kalman 

2.2  X  10“12 

1.3  X  10'6 

2.5  X  10'9 

5.3  X  10"5 

in  the  filter  models.  The  lack  of  process  noise  in  the  system  could 
result  in  artificially  low  estimates  of  the  covariance  matrix.  This 
also  causes  the  filters  to  follow  the  model  and  ignore  the  measurement 
updates,  resulting  in  slow  convergence. 

Also  noted  in  Table  II  is  the  slightly  better  performance  of  the 
Bayes  and  Kalman  filters  compared  to  the  Least  Squares  filter.  This  is 
in  part  due  to  the  basic  operation  of  the  filters.  The  Least  Squares 
filter  is  a  batch  estimator  and  updates  the  estimate  after  ten  measure¬ 
ment  updates.  The  filter  minimizes  the  residuals  over  the  last  ten 
measurements.  The  Bayes  and  Kalman  filters  are  recursive  filters  and 
minimize  the  residuals  of  the  last  measurement  only. 

Finally,  the  data  in  Table  II  indicates  that  the  Kalman  and  Bayes 
filters  provide  similar  results.  This  similarity  will  be  investigated 
further  in  the  following  sections. 

State  Estimate  Errors 

Another  measure  of  filter  performance  is  the  error  in  the  state 
estimates.  The  errors  evaluated  in  this  section  are  the  mean  error, 
the  standard  deviation,  and  the  root  mean  squared  (rms)  error. 


36 


[jjuiau  mpfj 


The  mean  error  given  in  the  previous  section  is 

1  N 

® j  “  jT  ^  (xT~xi ) i  (62) 

i=l  J 

This  provides  the  mean  error  of  the  jth  element  based  on  N  runs.  As 
stated  in  the  previous  section,  fifteen  runs  are  used  for  the  Monte 
Carlo  analysis. 

The  standard  deviation  is  the  square  root  of  the  variance  defined 
in  the  previous  section  and  the  rms  error,  a,  is  the  square  root  of  the 
mean  squared  error  defined  in  the  previous  section. 

Discussions  of  these  errors  for  each  satellite  case  are  presented 
in  the  following  sections. 

Case  I 

This  orbit  is  used  to  verify  filter  performance  under  normal 
operations.  A  measurement  update  is  input  every  10  time  units,  a  little 
more  than  two  hours  (133  min).  Twenty  measurement  updates  are  processed 
for  filter  evaluation. 

Plots  of  mean  error  and  mean  error  plus  and  minus  rms  error  are 
presented  for  comparison.  Plots  of  mean  error  plus  and  minus  standard 
deviation  may  be  more  meaningful;  however,  for  nearly  zero-mean  errors, 
there  is  little  difference  between  rms  error  and  standard  deviation. 

The  following  plots  show  that  the  errors  are  nearly  zero-mean. 

Figures  10-12  are  plots  of  the  mean  error  and  mean  error  plus  and 
minus  rms  error  of  the  semi -major  axis  for  the  three  filters.  The  mean 
errors  are  nearly  zero-mean.  The  Bayes  and  Kalman  Filters  demonstrate 
a  large  initial  transient  in  the  estimates.  The  mean  error  minus  the  rms 
error  is  nearly  800kmwhich  is  unacceptable  for  this  problem.  Both 
filters  recover  from  these  transients  to  errors  of  less  than  50  km. 
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Similar  results  are  observed  in  the  plots  for  mean  anomaly  for  the 
Bayes  and  Kalman  Filters  (Figs  13  and  14). 

The  large  initial  transients  could  be  due  to  two  possible  problems. 
First,  there  may  be  observability  problems  in  estimating  the  derivative 
of  the  semi -major  axis  a.  The  estimates  of  a  in  both  the  Bayes  and 
Kalman  Filters  tend  to  track  noise  in  the  data.  Second,  the  filters' 
initial  estimate  of  the  covariance  matrices  could  be  tuned  to  eliminate 
the  large  initial  transients. 

Another  difficulty  is  the  appearance  of  a  bias  in  the  estimation  of 
the  inclination  (Figs  15-17).  Also  the  Bayes  ad  Kalman  Filters  appear 
to  be  diverging  in  the  estimate  of  the  inclination.  These  problems 
will  be  investigated  in  Case  II  with  a  longer  time  scale. 

The  plots  of  mean  error  versus  time  for  the  remaining  elements 
are  contained  in  Appendix  F.  A  comparison  of  the  estimation  results  of 
the  Bayes  and  Kalman  Filters  shows  little  or  no  difference.  These 
similarities  are  also  demonstrated  in  Table  III,  which  presents  the 
maximum  mean  error,  maximum  rms  error,  and  maximum  standard  deviation 
for  each  element.  Due  to  their  similarities  subsequent  results  will 
be  presented  for  only  the  Bayes  or  Kalman  Filter,  but  not  both. 

Case  II 

This  orbit  is  input  to  evaluate  the  filter  performance  with  a 
reduced  rate  of  measurement  updates.  About  two  measurements  per  day 
are  input  to  the  filters.  Also,  the  semi -major  axis  is  two  earth  radii, 
so  atmosphere  effects  will  not  be  present.  Only  the  Least  Squares 
and  Bayes  Filters  are  evaluated  for  this  case. 

Figures  18  and  19,  plots  of  the  mean  error  in  the  semi-major  axis 
and  mean  anomaly  for  the  Least  Squares  Filter,  indicate  a  periodic 
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Kalman  Filter  -  Case 


Least  Squares  Filter  -  Case 


Fig  17.  ej  +  oj  vs  t,  Kalman  Filter  -  Case 


TABLE  III 


Maximum  Estimation  Errors  -  Case  I 


Least  Squares  Filter 


Element 

eft  (mean  error) 

a|vj  (rms  error) 

Max  Std  Dev 

a  (km) 

7.6 

25.0 

22.6  km 

e 

0.0009 

0.0013 

0.0012 

i  (deg) 

0.03 

0.03 

0.02 

oj  (deg) 

0.74 

0.81 

0.34 

£i  (deg) 

0.057 

0.069 

0.040 

M  (deg) 

1.26 

3.0 

2.90 

Bayes  Filter 

Element 

efl  (mean  error) 

om  (rms  error) 

Max  Std  Dev 

a  (km) 

286  (32)1 

486  (42)1 

395  (15)1 

e 

0.0008 

0.0028 

0.00265 

i  (deg) 

0.14 

0.24 

0.18 

u  (deg) 

1.1 

1.9 

1.5 

n  (deg) 

0.080 

0.11 

0.08 

M  (deg) 

8.6  (0.20)1 

17.2  (0.65)1 

'\ 

15.1  (0.63)1 

»  t 

Kalman 

Filter 

Element 

e^  (mean  error) 

cm  (rms  error) 

Max  Std  Dev 

a  (km) 

281  (36)1 

485  (42)1 

395  (15)1 

e 

0.0009 

0.0028 

0.00245 

i  (deg) 

0.13 

0.24 

0.18 

to  (deg) 

1.1 

1.9 

1.5 

n  (deg) 

0.08 

0.11 

0.10 

M  (deg) 

8.6  (0 . 1 7) 1 

17.2  (0.65)1 

15.1  (0.63)1 

Note  1: 

Maximum  values  following  initial  transient. 
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variation  with  a  period  of  about  500  TU,  about  five  days.  This  is  the 
frequency  of  the  estimate  update  for  this  filter.  The  Bayes  Filter 
estimates,  which  are  updated  every  50  TU,  about  twelve  hours,  do  not 
show  similar  trends  (Figs  20  and  21).  These  two  plots  again  display 
large  initial  transients  in  the  semi-major  axis  and  the  mean  anomaly. 
The  large  transients  are  in  the  rms  error  and  not  in  the  mean  error. 

Figure  22,  mean  error  in  inclination  versus  time  for  the  Bayes 
Filter,  shows  that  the  filter  is  again  diverging  in  the  estimate  of 
the  inclination. 

The  plots  of  mean  error  versus  time  for  the  remaining  elements 
are  contained  in  Appendix  F.  The  maximum  mean  errors,  rms  errors,  and 
standard  deviations  for  Case  II  are  shown  in  Table  IV. 

Case  III 

This  orbit  is  used  to  evaluate  filter  performance  with  large 
perturbations  due  to  air  drag.  The  Least  Squares  and  Kalman  Filters 
were  evaluated  for  this  orbit.  The  perigee  altitude  of  this  orbit  is 
200  km  and  the  apogee  altitude  is  1700  km.  The  drag  coefficient,  B,  of 
the  satellite  is  increased  from  0.01  to  0.1  Kg/m^.  The  non-circular 
orbit  is  used  to  provide  variations  in  the  air  drag  and  thus  the 
derivative  of  the  semi -major  axis.  Plots  of  mean  error  and  mean  error 
plus  and  minus  rms  error  in  the  semi -major  axis  verus  time  for  the 
two  filters  are  shown  in  Figures  23  and  24.  Similar  plots  for  the  mean 
anomaly  are  shown  in  Figures  25  and  26. 

Table  V  contains  the  maximum  mean  errors,  rms  errors,  and  standard 
deviations  for  each  element  for  both  filters. 


Fig  18.  ea  +  aa  vs  t,  Least  Squares  Filter  -  Case  II 


Fig  20.  ea  +  aa  vs  t,  Bayes  Filter  -  Case  II 
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TABLE  IV 


Maximum  Errors  -  Case  II 
Least  Squares  Filter 


Element 

e^  (mean  error) 

(rms  error) 

Max  Std  Dev 

a  (km) 

10.8 

40.8 

40.4 

e 

0.0014 

0.0017 

0.0010 

i  (deg) 

0.09 

0.032 

0.031 

w  (deg) 

0.26 

0.45 

0.36 

n  (deg) 

0.12 

0.13 

0.04 

M  (deg) 

1.09 

7.02 

7.02 

Bayes 

Filter 

Element 

e^  (mean  error) 

cm  (rms  error) 

Max  Std  Dev 

a  (km) 

30.0  (4.7)1 

212  (5.9)1 

204  (1.85)1 

e 

0.0006 

0.0036 

0 .0012 

i  (deg) 

0.08 

0.12 

0.09 

u  (deg) 

1.1 

1.27 

0.69 

«  (deg) 

0.014 

0.044 

0.043 

M  (deg) 

1.26  (0.23)1 

10.1  (0.50)1 

10.0  (0.38)1 

Note  1: 

Maximum  values  following 

initial  transients. 
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Least  Squares  Filter  -  Case  III 
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,  Least  Squares  Filter  -  Case  III 


Fig  26.  eM  +  cm  vs  t,  Kalman  Filter  -  Case  III 


TABLE  V 


Maximum  Errors  -  Case  III 


Least  Squares  Filter 


Element 

eM 

°M 

Max  Std  Dev 

a  (km) 

4.33 

16.0 

15.0 

e 

0.001 

0.001 

0.001 

i  (deg) 

0.003 

0.12 

0.03 

O)  (deg) 

0.74 

0.79 

0.33 

n  (deg) 

0.08 

0.09 

0.04 

M  (deg) 

0.86 

2.50 

2.46 

Kalman  Filter 

Element 

eM 

°M 

Max  Std  Dev 

a  (km) 

126  (12. 6)1 

588  (19. I)1 

580  (14. 0)1 

e 

0.0021 

0.0036 

0.0030 

i  (deg) 

0.20 

0.29 

0.21 

w  (deg) 

1.43 

1.98 

1.32 

fi  (deg) 

0.09 

0.14 

0.11 

M  (deg) 

5.27  (0.26)1 

36.7  (0 . 60) 1 

36.1  (0 . 57 ) 1 

Note  1: 

Maximum  values  following 

initial  transients. 
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Case  IV 

This  orbit  is  used  to  demonstrate  the  effects  of  singularities  on 
the  estimators  due  to  undefined  components  in  the  classical  elements. 

The  first  singularity  is  due  to  a  circular  orbit.  For  a  circular  orbit, 
the  argument  of  perigee  is  not  defined  since  the  perigee  itself  is  not 
defined.  Subsequently  the  mean  anomaly  is  not  defined.  A  zero  inclina¬ 
tion  orbit  presents  a  singularity  in  the  right  ascension  of  the  ascend¬ 
ing  node  since  the  ascending  node  is  not  defined. 

It  should  be  noted  that  these  singularities  are  actually  due  to  the 
transformation  from  position  and  velocity  vectors  to  classical  elements, 
and  are  not  singularities  in  the  orbit  itself.  The  orbit  is  still 
defined  by  six  elements;  but,  the  classical  elements  do  not  uniquely 
define  the  orbit. 

All  the  filters  have  difficulty  with  this  orbit.  The  major  problem 
is  estimating  the  mean  anomaly.  Since  this  quantity  is  not  uniquely 
defined,  large  residuals  occur.  The  filters  attempt  to  zero  out  these 
residuals  by  adjusting  the  mean  motion,  which  in  turn  causes  the  semi¬ 
major  axis  to  vary.  The  estimate  of  the  semi -major  axis  diverges  rapidly; 
and,  in  the  case  of  this  orbit,  the  filters  eventually  estimate  a  negative 
semi -major  axis,  halting  program  execution. 

This  represents  a  major  problem  for  the  filters  used  in  this  study. 

A  possible  solution  to  this  problem  is  to  use  equinoctal  orbital  elements 
as  the  filter  states  [Ref  9;  44].  These  elements  do  not  have  singular¬ 
ities  due  to  circular  orbits  or  zero  inclination  orbits.  Classical 
elements  could  still  be  used  as  measurement  updates.  The  only  change 
would  be  the  observation  relationship,  H.  The  actual  estimator  equations 
would  not  be  changed. 
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Case  V 


This  orbit  presents  observation  problems  to  the  estimators.  Due 
to  the  highly  eccentric  orbit,  the  satellite  is  beyond  the  range  of 
the  radars  most  of  the  time;  therefore,  large  data  gaps  are  present  in 
the  data. 

Attempting  an  observation  every  two  hours,  only  one  observation 
per  day  is  obtained  due  to  the  range  restrictions.  Due  to  the  long 
running  time  of  the  truth  model  to  provide  data  for  the  filters,  a 
Monte  Carlo  analysis  for  this  case  was  not  completed.  A  single  run 
using  the  Kalman  and  Least  Squares  Filters  shows  no  difficulty  in 
estimating  the  orbital  elements  for  this  highly  eccentric,  low  perigee 
orbit.  Table  VI  shows  the  estimation  errors  for  the  two  filters  after 
twenty  measurement  updates,  about  twenty  days. 


TABLE  VI 

Estimation  Errors  for  Case  V 


Element 

Least  Square 

Kalman 

a 

25.5  km 

57.0  km 

e 

0.001 

0.0008 

1 

0.006° 

0.03° 

u 

0.20° 

0.12° 

fl 

0.97° 

0.02° 

M 

3.15° 

2.75° 

It  should  be  noted  that  these  errors  are  only  for  one  run  and 
are  not  the  result  of  a  Monte  Carlo  analysis. 
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Discussion 

As  seen  in  the  previous  examples,  all  filters  encounter  problems 
in  estimating  the  elements  and  the  errors. 

First,  estimation  was  not  possible  with  singularities  in  the 
classical  elements.  Second,  the  Bayes  and  Kalman  Filters  demonstrate 
large  initial  transients  in  the  estimation  problem.  These  initial 
transients  are  reduced.  These  transients  are  probably  due  to  observa¬ 
bility  problems  in  the  derivative  of  the  semi-major  axis  and  an 
improperly  tuned  filter. 

Next,  the  in-track  errors  in  the  estimates  are  too  large  for  SDC 
application.  This  may  be  due  to  using  too  large  errors  in  the  initial 
radar  data  or  failure  of  the  filters  to  converge  within  the  specified 
number  of  iteration. 

Finally,  the  Bayes  and  Kalman  Filters  track  noise  in  the  estimates 
of  the  derivatives  of  the  semi-major  axis  and  the  eccentricity.  The 
wrong  sign  on  these  values  would  greatly  effect  the  ability  of  the 
filters  to  propagate  the  states  forward  and  predict  the  satellite's 
position  for  radar  acquisition. 

In  spite  of  the  problems  listed  above,  the  potential  of  using 
two-body  elements  has  been  demonstrated  in  this  study.  If  the  prob¬ 
lems  listed  above  can  be  eliminated,  any  of  the  filters  evaluated 
could  be  used  by  SDC  for  orbital  estimation. 

Another  characteristic  noted  in  the  filter  performance  is  the 
tendency  of  the  gain  matrices  in  the  Kalman  and  Bayes  Filters  to 
approach  steady  state  for  a  particular  orbit. 


As  shown  in  Chapter  III,  the  gains  for  the  Bayes  and  Kalman  are 
re-evaluated  for  every  measurement  update.  After  several  measurements 
the  gain  matrices  change  very  little  from  update  to  update.  This  is 
due  to  the  nearly  constant  values  of  the  orbital  elements.  If  a  pre¬ 
determined  gain  can  be  stored  in  the  computer  a  significant  reduction 
in  computer  time  can  be  achieved.  In  the  case  of  the  Kalman  Filter,  no 
matrix  inversions  would  be  required.  The  covariance  update  would  be 
given  by 

P(+)  =  P(-)  -  KH*P(-)  (63) 

and  the  change  in  the  state  estimate  would  be  given  by 

6x(+)  =  6x(l-l)+  K(r.  -  H$6x(-))  (64) 

The  calculation  of  the  gain 

K  =  P(-)$T  HT(Q  +  H$P  H1*)'1  (65) 

would  not  be  required  for  every  measurement  update. 
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VI.  Conclusions  and  Recommendations 


Conclusions 

The  following  conclusions  are  drawn  from  the  filter  performance 
analysis  completed  in  this  study: 

1.  The  feasibility  of  using  two-body  classical  orbital  elements 
as  measurement  updates  for  orbital  element  estimation  is  demonstrated 
by  all  three  filters. 

2.  All  three  filters  diverge  when  estimating  the  elements  of  a 
circular  orbit. 

3.  In-track  errors  in  the  filter  estimates  are  too  large  for  SDC 
requirements  and  techniques  must  be  developed  to  reduce  the  error. 

4.  All  filters  underestimate  the  covariance  of  the  estimates. 

5.  The  Kalman  and  Bayes  Filters  are  similar  and  differences  in 
their  performance  can  not  be  detected. 

6.  All  three  filters  track  noise  in  the  estimation  of  the 
derivative  of  the  semi -major  axis. 

7.  The  Bayes  and  Kalman  Filters  have  large  initial  transients 
in  the  estimation  of  the  semi -major  axis  and  the  mean  anomaly.  These 
transients  are  reduced  with  further  updates. 

Recommendations 

1.  Change  the  filter  states  to  equinoctal  orbital  elements  to 
eliminate  problems  due  to  singularities  in  the  classical  elements. 

2.  Add  process  noise  to  the  filter  dynamics  and  tune  the  filters 
to  improve  covariance  estimates  and  remove  initial  transients. 
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3.  Re-program  the  Least  Squares  Filter  to  improve  computational 
efficiency  by  taking  advantage  of  symmetry  and  zeroes  in  the  matrices. 

4.  Modify  the  Bayes  and  Kalman  Filters  to  provide  batch  estimates 
of  the  derivative  of  the  semi -major  axis. 

5.  Obtain  actual  satellite  data  for  further  filter  evaluation. 

6.  Investigate  the  possibility  of  using  a  constant  gain  matrix 
for  each  orbit  in  the  Bayes  and  Kalman  Filters. 

7.  Investigate  other  methods  of  estimating  the  air  drag  such  as 
parameter  identification. 

8.  Change  the  Bayes  and  Kalman  Filters  from  infinite  memory 
filters  to  finite  or  fading-memory  filters. 

9.  Investigate  the  estimation  of  the  two-body  orbital  elements 
from  radar  data  to  gain  insight  into  evaluating  the  measurement  noise 
used  in  the  filters. 
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Appendix  A 

Coordinate  Frame  Transformations 


There  are  two  coordinate  frame  transformations  required  in  this 
study.  The  first  is  a  transformation  from  the  perifocal  to  the  geocentric- 
equatorial  coordinate  frame  [Fig  27].  It  is  used  in  the  generation  of  the 
satellite's  initial  position  and  velocity  vector  from  the  input  orbital 
elements.  The  transformation  is  given  by  [Ref  3,  80-83]: 


r  - 
X 

Xu 

Y 

ii 

1  1 

30 

1 _ 1 

Yco 

Z 

Zu 

where 

R  is  the  direction  cosine  matrix  between  the  two  frames. 

The  elements  of  R  are 

R11  =  cos  n  cos  a)  -  sin  fi  sin  u  cos  i 

R]2  =  -cos  n  sin  w  -  sin  n  cos  w  cos  i 
R]3  =  sin  n  sin  i 

Rgi  =  sin  n  cos  00  +  cos  n  sin  w  cos  i 

R22  *  -sin  o  sin  (*)  +  cos  n  cos  u  cos  i 

R23  =  -cos  Q  sin  i 

R3]  =  sin  u>  sin  i 

R32  =  cos  jo  sin  i 

R33  =  cos  i 

The  angles  i,  o  ,  w  are  the  inclination,  right  ascension  and 
argument  of  perigee. 
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The  second  coordinate  frame  transformation  is  from  the  local 
horizon  coordinate  frame.  This  transformation  is  used  to  rotate  the 
geopotential  acceleration  terms  into  the  geocentric-equatorial  frame 
for  integration.  The  rotation  is  given  by  [Ref  9;  91]. 


where 


[T] 


where 


X 
Y 
Z 

COS  a  COS  6 

sin  a  cos  6 
sin  6 


-sin  a 

COS  a 
0 


a  is  the  geocentric  right  ascension 
6  is  the  geocentric  declination 


(63) 

-COS  a  sin  <5 
-sin  a  sin  6 
cos  6 
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Appendix  B 

Position  and  Velocity  Data  from  Classical  Elements 


In  this  study,  a  satellite's  initial  position  and  velocity  vectors 
are  derived  from  an  input  orbit.  The  equations  needed  for  this  trans¬ 
formation  are  presented  below  [Ref  3,  72-73] : 


P  =  a(l-e2) 

(65) 

E  -  e  sin  E  =  M 

(66) 

_  e  -  cos  E 
cos  v  e  cos  E  -  1 

(67) 

r  =  p 

1  +  e  cos  v 

(68) 

sin  v  =  a*^ “e2  sin  E 
r 

(69) 

A  A 

r  =  r  cos  v  P  +  r  sin  Q 

(70) 

v  [-sin  v  P  +  (e  +  cos  v  )  Q] 

n 

(71) 

The  position  and  velocity  vector  are  then  transformed  from  the 
perifocal  to  the  geocentric-equatorial  coordinate  frame. 

It  should  be  noted  that  equation  (66)  is  known  as  Kepler's  equation 
and  a  closed  form  solution  for  the  eccentric  anomaly,  E,  cannot  be 
obtained.  A  Newton  iteration  scheme  is  used  to  evaluate  E,  using 
Eo  =  M  as  an  initial  guess. 
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Appendix  C 

Geopotential  Constants  and  Legendre  Polynomials 

Listed  below  are  the  geopotential  coefficients  for  the  zonal  and 

sectoral  harmonies  used  in  this  study  [Ref  6]: 

J2  =  1.082637  X  10'6 

(72) 

J3  =  -2.541  X  10'6 

(73) 

J4  =  -1.618  X  10‘6 

(74) 

C22  =  1-5362  X  10"6 

(75) 

S22  =  -.86358  X  10“6 

(76) 

Also,  the 

following  Legendre  polynomials  are  required  for  evalua 

tion  of  the  perturbing  accelerations  due  to  the  geopotential 

[Ref  9; 

90-91]: 

P0  =  i 

(77) 

Pi  =  sin  6  i 

(78) 

* 

1 

P2  =  ^C-l  +  3  sin2fi] 

(79) 

P3  =  -^[-2  sin  6  +  5  P2  sin  6] 

(80) 

p4  -  ^-3  P2  +  7  P3  sin  6] 

(81) 

p|  =  3  cos2  <5 

(82) 

cos  <5 

=  3  sin  6  cos  <5 

(83) 

cos  6  Pg1 

-  sin  6  cos  6  +  3  cos  6  P2 

(84) 

COS  6  P^' 

=  sin  6  cos  6  Pg1  +  4  cos  6  P3 

(85) 

COS  6  P%' 

p! 

’  -2  s,n  5  cos  s 

(86) 

(86) 


where 


6  is  the  geodetic  declination  of  the  satellite. 
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Appendix  D 

Orbital  Elements  from  Position  and  Velocity  Data 

The  estimators  in  this  study  use  classical  orbital  elements  as 
measurement  updates.  These  elements  are  calculated  from  the  position 
and  velocity  data  from  the  truth  model.  Prior  to  calculating  the 
elements,  noise  is  added  to  the  data.  The  following  equations  are  used 
to  obtain  classical  elements  from  thegeocentric-equatorial  position  and 
velocity  data  [Ref  3;  62-63]: 


h  = 

L  x  v. 

(87) 

n  = 

k  x  h 

(88) 

e  = 

J[(v-f)r.(p.y)v] 

(89) 

e  = 

lei 

(90) 

P  - 

h2/p 

(91) 

P 

(92) 

a 

TT^eZy  , 

cos 

i  =  T  0  <  1  <  180 

(93) 

cos 

Q  =  — -  If  n-j  >  0,  n  <  180° 
n  J 

(94) 

cos 

a)  =  -  ’  -  If  ek  >  0,  a»  <  180° 
ne  * 

(95) 

e  •  r 

(96) 

cos 

v  =  -  If  r  •  v  >  0,  v0  <  180° 

cos 

r  _  e  +  cos  v 

1  +  e  cos  v 

(97) 

M  =  E  -  e  sin  E 

(98) 
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where 

is  the  angular  momentum  vector 
n^  is  the  nodal  vector 

The  eccentric  anomaly,  E,  should  be  in  the  same  half-plane  as  the 
true  anomaly,  v0. 

This  method  provides  a  set  of  osculating  elements  which  match  the 
orbit  at  the  time  the  position  and  velocity  is  generated. 


Appendix  E 

Derivation  of  State  Transition  Matrix 


The  state  transition  matrix,  $  ,  is  calculated  from  equation  (16) 
in  Chapter  III. 


or. 


9  = 


where 


3fl/axi  3X2 
af2/axi  *  '  • 


3fl_ 

3X8 

3f2 

3*8 


L 


3f8/axi 


3fg 

3*8 


(99) 


x(t)  =  f(x(t0),  t) 

is  an  approximate  solution  to  the  non-linear  system. 

The  equations  for  the  state  dynamics  presented  in  Chapter  III 


(equation  (16))  are  listed  below: 

f]  =  a(t)  =  a0  +  a(t  -  t0)  (100) 
f2  =  e(t)  =  e0  +  e(t  -  t0)  (101) 
f3  =  1(t)  =  io  (102) 
U  =  «(t)  =  o)0  +  ;(t  -  t0)  (103) 
f5  =  0(t)  =  00  +  „(t  -  t0)  (104) 
f6  =  M(t)  =  M0  +  n|c(t  -  t0)  +  j  n(t  -  t0)2  (105) 
fj  =  a(t)  =  a0  (106) 
f8  =  e(t)  =  et>  (107) 


where 


f 
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=  n  +  e 


(108) 


-3n  J?  re2  R 

a  =  - — •  (f  sin2  i-2) 

2a2(l-e2)2  '2 

(109) 

-3n  J2  re2 

a  =  - — -  cos  i 

2a2(i-e2)2 

(no) 

n  =  ^  H 
n  2  a  a 

(111) 

-3n  J2re2  ,3  .  ?  .  -.x 
e  =  .  ,/9  (0  sin‘  1-1) 

2a2(i-e2)3/2  2 

(112) 

Evaluation  of  the  partial  derivative  provides  the  following 
elements  to  the  <f  matrix: 


$ii  =  1  1=1,8 

(113) 

*1,7  *  *2,8  =  t  -  t0 

(114) 

..  .  _  21  n  J2  re2  ,5  4  ow*  *  x 

‘  4  a3(lHS2)2  2  s,n  *  *»> 

(115) 

-6n  J9  r  2e  c 

^■a2(,.e2,3  <f «M*  -  *•> 

(116) 

*4.3  =  2a2(l-e2)2  (S’"  ’  C°S  '  ‘°) 

(117) 

21 n  J2  rp2  .  ,  ,  . 

*5'’  =  4a2(,-e2)2  (C°S  *  (t  -  ‘o’ 

(118) 

6n  Jo  rp2e 

(119) 

*5,2  =  9  0  (cos  i  (t  -  t0) 

a*(l-e^)J 

3n  J2  re2  ...  v 

*5-3  '  2a2(l-e2)2  S’"  ’  lt  "  t0> 

(120) 

^6,1  ‘  2  a  (t  ‘  ^  +  4a3(1.e2)2  ^  S1n  1  1)(t  to) 

(121) 
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022) 

-9n  Jo  re2 

(123) 

*6,3  '  2a2(l-62)3«  S,n  '  C°S  '  “  '  to) 

♦6,7  "  Ta  (t  -  ‘o)2 

(124) 

The  remaining  terms  in  the  $  matrix  are  all  zero. 
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Appendix  F 

Plots  for  Filter  Performance  Analysis 

The  following  plots  are  presented  to  supplement  the  discussion 
and  data  presented  in  Chapter  V.  The  plots  are  the  mean  error  and 
mean  error  plus  and  minus  rms  error  versus  time.  As  stated  in 
Chapter  V,  the  rms  error  is  nearly  equal  to  the  standard  deviation. 


Fig  28.  Og  vs  N,  Least  Squares  Filter 


Fig  29.  c/  vs  N,  Least  Squares  Filter 
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vs  N 


Fig  35.  Op  vs  N,  Kalman  Filter 


Fig  36.  of  vs  N,  Kalman  Filter 
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Fig  37.  o,  vs  N,  Kalman  Filter 
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Fig  38.  o„  vs  N,  Kalman  Filter 


Fig  39.  ee  +  cre  vs  t,  Least  Squares  Filter  -  Case 


Least  Squares  Filter  -  Case 


Least  Squares  Filter  -  Case 
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Fig  44.  eu  +  a  vs  t,  Bayes  Filter  -  Case 


Bayes  Filter  -  Case 


Fig  46.  ee  +  oe  vs  t,  Kalman  Filter  -  Case 


Fig  48.  en  +  ou  vs  t,  Kalman  Filter  -  Case 


Fig  49.  ep  +  ae  vs  t.  Least  Squares  Filter  -  Case  II 


wo 
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Least  Squares  Filter  -  Case  II 


Fig  53.  ee  +  oe  vs  t,  Bayes  Filter,  Case  II 
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able  in  this  study.  All  three  filters  performed  adequately  with  the  following 
exceptions.  The  filters  diverge  when  singularities  are  present  in  the  classi¬ 
cal  orbital  elements.  All  three  filters  underestimate  the  covariance  of  the 
estimates.  Finally,  the  filters  track  noise  in  the  estimate  of  the  derivative 
of  the  semi -major  axis,  degrading  the  prediction  capabilities  of  the  filters. 
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